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We study quantum turbulence in trapped Bose-Einstein condensates by numerically solving the 
Gross-Pitaevskii equation. Combining rotations around two axes, we successfully induce quantum 
turbulent state in which quantized vortices are not crystallized but tangled. The obtained spectrum 
of the incompressible kinetic energy is consistent with the Kolmogorov law, the most important 
statistical law in turbulence. 



PACS numbers: 03.75.Kk 05.30.Jp 47.32.C- 47.37.+q 

The study of turbulence has a very long history, go- 
ing back at least to Leonardo da Vinci, and understand- 
ing and controlling turbulence are great dreams of sci- 
ence and technology. Classical turbulence (CT) exhibits 
highly complicated configurations of eddies. Many stud- 
ies have been devoted to the dynamical and statistical 
properties of CT after Kolmogorov's pioneering work 
[l|, l2| on flow at very high Reynolds number, namely, 
fully developed CT. The characteristic behavior of CT 
has been believed to be sustained by the Richardson cas- 
cade of eddies from large to small scales. However, in 
CT, there is no universal way to identify each eddy, be- 
cause they continue to nucleate, diffuse, and disappear. 
As a result, many aspects of CT are still not perfectly 
understood. 

Turbulence is also possible in superfluids, such as the 
superfluid phases of 4 He and 3 He. Such quantum turbu- 
lence (QT) consists of definite topological defects known 
as quantized vortices and has recently attracted interest 
as a way to better understand turbulence [3| . 

Superfluid 4 He has been extensively studied, in par- 
ticular with relation to quantized vortices Below 
the lambda temperature T\ = 2.17 K, liquid 4 He en- 
ters the superfluid state through Bose-Einstein conden- 
sation. The hydrodynamics of superfluid 4 He is strongly 
influenced by quantum effects; any rotational motion is 
sustained by quantized vortices with quantized circula- 
tion k = h/m, where m is the particle mass. There are 
two typical cooperative phenomena of quantized vortices. 
One is a vortex lattice under rotation in which straight 
quantized vortices form a triangular lattice along the ro- 
tation axis p|. The other is a vortex tangle in QT in 
which vortices become tangled in a flow [1,13|- 

QT has been studied as a problem in low temperature 
physics since its discovery some 50 years ago. Its study 
has recently entered a new stage beyond low temperature 
physics. One of the main motivations of recent studies is 
to investigate the relationship between QT and CT. Some 
similarities between the two types of turbulence have 
been experimentally observed in superfluid 4 He and 
3 He [HI [n|, and have been theoretically confirmed by 
numerical simulations of the quantized vortex-filament 
model fl2| and a model using the Gross-Pitaevskii (GP) 
equation [13L [13. Il5l Il6j. In particular, we have success- 
fully obtained the Kolmogorov law for QT, which is one 



of the most important statistical laws in CT |I7I| by a 
numerical simulation of the GP equation [3, Il5|. 

The similarity between QT and CT means that QT is 
an ideal prototype to study the statistics and vortex dy- 
namics of turbulence, because QT exhibits a real cascade 
process of quantized vortices. However, in superfluid he- 
lium, it is very difficult to experimentally control the tur- 
bulent state and determine the vortex configuration. 

Another important example of quantized vortices is 
magnetically or optically trapped atomic Bose-Einstein 
condensates (BECs) [l8j]. The characteristics of trapped 
BECs are as follows: (i) a BEC system is weakly inter- 
acting and can be easily treated theoretically, (ii) many 
physical parameters of BECs are experimentally control- 
lable, and (iii) various physical quantities such as the den- 
sity and phase of BECs can be directly observed, which 
is in stark contrast to superfluid helium systems. Quan- 
tized vortices can be considered to be holes of density 
and singularities of phase. Shortly after trapped BECs 
were first realized, experimental groups [HI, [2(| reported 
vortex lattice structures, as well as the crystallization 
dynamics of these structures under rotation. These dy- 
namics have been successfully confirmed quantitati vely 
by the numerical simulation of the GP equation [2l], |22| . 

However, in the experimental research of trapped 
BECs, another important phenomenon of quantized vor- 
tices, namely QT, has not been adequately studied. Not- 
ing that quantized vortices are observable and that al- 
most all physical parameters of trapped BECs are con- 
trollable, such systems are an ideal prototype for truly 
controllable QT, which is not possible for superfluid he- 
lium. QT in trapped BECs can be used to determine 
several details of the system, such as the distribution 
of vortex length, details on the cascade of vortices, the 
isotropy or anisotropy of vortex configuration and details 
on correlations among vortices related to eddy viscosity, 
as already considered for CT [TjJ . Clarifying any of these 
will allow the transition to QT to be considered a univer- 
sal phase transition. Therefore, research into QT offers 
the promise of greater advances in understanding turbu- 
lence than has been possible in past studies of turbulence. 

There are only a few theoretical pioneering works for 
QT in trapped BECs in which QT is realized in relax- 
ation processes from non-equilibrium initial state across 
the BEC critical temperature [23[ , or from vortex- free to 
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vortex lattice state under rotation [l6[. In this paper, 
we present a numerical simulation of the 3-dimensional 
GP equation. First, we present the more experimentally 
realizable steady QT in a trapped BEC by combining 
rotations around two axes. Second, we show that the 
spectrum of the incompressible kinetic energy (k) per 
unit mass obeys the Kolmogorov law 



1st rotation 



E'l in (k) = Ce 2/3 k- 5 / 3 



(1) 



Here, the energy spectrum is defined as E^ in = 
J dfci?J. in (fc), where k is the wave number from the Fourier 
transformation of the incompressible velocity field and e 
is the energy transportation rate from small to large k 
of the incompressible kinetic energy per unit mass. The 
Kolmogorov constant C is a dimensionless parameter of 
the order of unity in CT. We further obtain the length 
distribution of vortices with a scaling structure that re- 
flects the self-similarity in the Richardson cascade. 

In considering a trapped BEC system, we start from 
the GP equation 



[i- 7 (#|^,t) = 



2m 



+g\$>{x,t)\ 2 + U(x) - H(t) • L(x) 



$(*,*)• (2) 



Here $(x,t) = f{x,t)e 1 ^ x ' t ^ is the macroscopic wave 
function of the BEC, m is the particle mass, p is the 
chemical potential, L(x) = —ifvx x V is the angular mo- 
mentum, and g = Airh 2 a/m is the coupling constant with 
s-wave scattering length a. The trapping potential U (x) 
is given by a weakly elliptical harmonic potential: 



.[(i-<yi)(i-<y a )a^ 



U(x) 



+ (l + 6 1 )(l-5 2 )y 2 + (l + 6 2 )z% 



(3) 



where us is the frequency of the harmonic trap and the 
parameters Si and 5 2 exhibit elliptical deformation in 
the xy- and za;-planes. To develop the BEC to a tur- 
bulent state rather than a vortex lattice state, we com- 
bine two rotations along the z-and shown in 
Fig. [TJ The rotation vector fi(t) is given by fi(£) = 
(f2 x , O z sinO^i, J7 2 cosil^i), where Q z and £l x are the fre- 
quencies of the first and second rotation, respectively. 
There are two main advantages of using the combined 
rotation to study turbulence. First, we can do direct 
numerical simulations and directly compare them with 
experiments without any ambiguity. Second, it is pos- 
sible to freely control the state from non-turbulent vor- 
tex lattice to fully-developed turbulence by changing the 
ratio Vl x /Vl z . In the classical fluid system, Goto et al. 
have already adopted this combined rotation to study 
turbulence by using water in spinning sphere on a rotat- 
ing turntable [24j and reported the transition from rigid 
body rotation to non periodic turbulent motion of water. 
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FIG. 1: (color online). Schematic sketch of the rotation. The 
first rotation is applied along the z-axis and the second rota- 
tion is applied along the x-axis. 



The condensate density p(x,t) and the superfluid ve- 
locity v(x,t) are given by p(x, t) — f 2 (x, t) and v{x 1 t) = 
H/rriV<f>(x,t). The vorticity rotv(x,t) vanishes every- 
where in a single-connected region of the fluid; any ro- 
tational flow is carried only by quantized vortices in the 
core, of which $>(x,t) vanishes so that the circulation 
is quantized by n = 2wh/m. The vortex core size is 
given by the healing length £ = h/^/2mgp. In trapped 
BECs, the healing length depends on the position, be- 
cause the system is not uniform. In this work, we define 
the characteristic healing length £ = %/ ^/2mgpa with the 
condensate density po = p(x = 0) at the trap center. 

We take a system at very low temperatures. The phe- 
nomenological damping term j(x) simulates the effect of 
thermal excitations and is effective only at scales smaller 
than the core size of quantized vortices, as shown in our 
previous work on the numerical simulation of a coupled 
system using the GP equation and the Bogoliubov-de 
Gennes equation (see Fig. 2(b) in the reference [25j] ) . 
In this work, we adopt the Fourier transformed damping 
term j(k) = j(k) at the temperature T — 0.01T C given in 
our previous work, where T c is the BEC critical temper- 
ature of ideal Bose gas. Introduction of j(x) conserves 
neither the energy nor the number of particles. To avoid 
the inconservation of the number of particles, we consider 
the time dependence of the chemical potential so that the 
total number N — J dx |<i>(a:,£)| 2 can be conserved. 

To solve the GP equation ([2]) numerically with high 
accuracy, we use the pseudo-spectral method, applying 
the Chebyshev tau method in space with a Dirichlet 
boundary condition in a box (2o| containing 512 3 grid 
points. For the numerical parameters, we use the fol- 
lowing, taken from experiments on 87 Rb atoms (l8l. Il9j|: 
rn = 1.46 x 1CT 25 kg, a = 5.61 nm, N = 2.50 x 10 5 , 
us = 150 x 2ir Hz. The total volume of the numerical box 
is set to V = 14. 3 ^m 3 . 

We start from a stationary solution without rotation 
and elliptical deformation. At t = 0, we turn on the 
rotation fl x = fl z = 0.6 and elliptical deformation 
Si = 5 2 = 0.025, and numerically calculate the time 
development of the GP equation @ using the 4th or- 
dered Runge-Kutta method [26J] with a time resolution 
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of At = 1 x 10~ 4 w. 

We calculate the total compressible and incompressible 
kinetic energy per unit mass E^ in (t) and E^ in (t) denned 
by 



1 

27V 



dx{[p(x,t)rr 



(4) 



Here, p(x,t) = h/mf(x,t)'V4 > (x,i), [•■•] c denotes the 
compressible part V x [• • • ] c — and [•••]' denotes the 
incompressible part V • [• • • J 1 = 0. We also investigate the 
anisotropy of the system by defining a parameter a(t) as 



A(t) = J dx |V$(a?,t)| 2 , 



m = E 



A(t) 



a(t) = 



A{ty 



(5a) 
(5b) 
(5c) 



for xi 



2,3 



{x,y, z}. Figures HJa) and (b) show the 



time development of E?£. 1 (t) and a(t) respectively. After 
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FIG. 2: Time development of E^(t) (a) and a(t) (b) 

tui ~ 150, a[t) becomes small showing that the BEC 
recovers isotropy, and E^ in (t) becomes almost steady, 
which means that isotropic steady turbulence is realized 
at tu; > 150. The steady turbulence is sustained by the 
balance between the large-scale energy injection due to 
the rotation and the small-scale dissipation. The time of 
tu> = 150 corresponds to t = 0.1 sec which is sufficiently 
shorter than the actual lifetime of trapped BECs in ex- 
periment. E^ in (t) is always much larger than E^ in (t) and 
dynamics of the BEC is dominated by vortices rather 
than compressible excitations in the BEC. 

To confirm that the system is genuine turbulence, we 
calculate the spectrum El in (k,t) and the flux U(k,t) 
from small to large k of the incompressible kinetic en- 
ergy El in (t). II(fc, t) can be obtained by considering the 
scale-by-scale energy budget equation for the GP equa- 
tion and given as (see Eq. (35) in the reference (l5l |) 



n(fc,t) = 



dx L k [{p(x,t) ■ Wv(x,t)Y] 
■ L k [{p(x,t)Y]. (6) 



Here, L k is the operator for the low-pass filter: 

£*[*(*)] = i £ Jdx'e ik < x - x 'Ks(x), (7) 
|fe|<fc 

for an arbitrary function s(x). Figure [3] (a) shows 
E^ in (k,t) and H(k,t) for the turbulent state. E^ in (k,t) 
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FIG. 3: (color online), (a) Wave number dependence of 
El. in (k,t) and H(k,t). (b) Vortex Length distribution n(l)Al 
inside the Thomas-Fermi radius Rtf- El. in (k,t), II(fc,t) and 
n(l)Al are obtained from an ensemble average of 25 turbu- 
lent states (ta> > 300). In both figures, a tra p = y/h/mu) is the 
characteristic scale of the trap. 

in this QT satisfies the Kolmogorov law of Eq. (JXJ) in 
the inertial range 2it/Rtf < k < 27r/£, where Rtf — 
y/2(j,(t = 0)/mw 2 is the Thomas-Fermi radius and repre- 
sents the largest scale in the BEC [27j . Furthermore, the 
energy flux is nearly constant value U(k,t) ~ \Ahia 2 jm 
in the inertial range, supporting that the incompress- 
ible kinetic energy steadily flow in wave number space 
through the Richardson cascade at the constant energy 
transportation rate e = H(k,t) in Eq. ([1]). Using this 
e, we obtain the Kolmogorov constant C = 0.25 ± 0.2 
which is smaller than that in CT and consistent with our 
previous work for QT in the uniform system [Til . ll5T ] . 

To investigate the relation between the Kolmogorov 
law and the Richardson cascade, we calculate the vortex 
length distribution n(l)Al inside the condensate (Fig. [3] 
(b)), where n(l)Al represents the number of vortices with 
length from Z to Z + Al. At the turbulent state, n(l)Al 
obeys the scaling property n(l)Al oc Z~" for 2ir£ < I < 
2ttR-jf . This reflects the self-similar Richardson cascade 
in which large vortices entering the condensate from the 
surface [2l], [23] are divided into smaller vortices, which 
is first confirmed in turbulence with the framework of 
the Gross-Pitaevskii equation. The scaling exponent a 
is close to unity, which is consistent with those given by 
Araki et al. (a ~ 1.34) [l2| and Mitani et al. (a ~ 1) 

To visualize the turbulence, we plot the isosurface of 
the condensate density p(x,t) and the spatial distribu- 
tion of the vortices inside the condensate in Fig. [4] (a)- 
(f). At tu> = 10, the surface of the BEC becomes unstable 
(Fig. 2] (a) and (d)), and vortices appear in the BEC at 
tu = 50 (Fig. H](b) and (e)). FiguresH (c) and (f) shows 
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FIG. 4: (color online). Isosurface plots of 5% of the maximum 
condensate density (a)-(c) and configuration of quantized vor- 
tices inside the Thomas- Fermi radius _Rtf (d)-(f). (a), (d) 
tu = 10; (b), (e) tu = 50; (c), (f) tu = 300. The method for 
identifying vortices in (d)-(f) is the same as that in Fig. 7 in 
the reference [l5j]. 

QT with no crystallization but with highly tangled quan- 
tized vortices at tu> = 300. 

In conclusion, using a numerical simulation of the GP 
equation, we have induced QT in a trapped BEC by com- 
bining rotations around two axes. The quantized vortices 



in the trapped BEC are not crystallized but tangled. In 
the inertial range, the spectrum of the incompressible ki- 
netic energy obeys the Kolmogorov law and the energy 
flux becomes constant value. We further obtained a vor- 
tex length distribution with a scaling property, support- 
ing the self-similar Richardson cascade of quantized vor- 
tices. The incompressible kinetic energy and its spectrum 
can be experimentally observed by measuring the density 
and phase of the BEC, according to Eq. ([4]). We antic- 
ipate the experimental realisation of QT in a trapped 
BEC and further advancement of the understanding of 
turbulence. 

The obtained energy spectrum in Fig. [3] (a) is not 
so destructive straight line and its consistency with 
the Kolmogorov law is incomplete. This inconsistency 
comes from the anisotropy of turbulence around the y- 
axis around which there is no rotation, and is resolved 
by other simulations of quantum turbulence of trapped 
BECs under the combined rotations around three axes. 
We will report on this study in near future. 
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